An elevational gradient in floral traits and pollinator assemblages in the Neotropical species Costus guanaiensis var. tarmicus in Peru

Abstract Different populations of plant species can adapt to their local pollinators and diverge in floral traits accordingly. Floral traits are subject to pollinator‐driven natural selection to enhance plant reproductive success. Studies on temperate plant systems have shown pollinator‐driven selection results in floral trait variation along elevational gradients, but studies in tropical systems are lacking. We analyzed floral traits and pollinator assemblages in the Neotropical bee‐pollinated taxon Costus guanaiensis var. tarmicus across four sites along a steep elevational gradient in Peru. We found variations in floral traits of size, color, and reward, and in the pollinator assemblage along the elevational gradient. We examined our results considering two hypotheses, (1) local adaptation to different bee assemblages, and (2) the early stages of an evolutionary shift to a new pollinator functional group (hummingbirds). We found some evidence consistent with the adaptation of C. guanaiensis var. tarmicus to the local bee fauna along the studied elevational gradient. Corolla width across sites was associated with bee thorax width of the local most frequent pollinator. However, we could not rule out the possibility of the beginning of a bee‐to‐hummingbird pollination shift in the highest‐studied site. Our study is one of the few geographic‐scale analyses of floral trait and pollinator assemblage variation in tropical plant species. Our results broaden our understanding of plant‐pollinator interactions beyond temperate systems by showing substantial intraspecific divergence in both floral traits and pollinator assemblages across geographic space in a tropical plant species.

Selection often changes floral traits in a way that enhances plant reproduction by increasing attraction and reward traits for pollinators, resulting in higher visitation rates, and by shaping traits in a way that improves the pollinator-flower fit, resulting in higher pollination efficiency (Stebbins, 1970). Thus, different plant populations of the same species adapt to local pollinators and diverge in floral traits accordingly (Anderson et al., 2010;Anderson & Johnson, 2008Galen, 1996;Johnson & Steiner, 1997;Maad et al., 2013;Medel et al., 2007;Nattero et al., 2011;Newman et al., 2015;Thompson, 2005). Understanding floral adaptation to pollinators can provide insight into major evolutionary processes involved in angiosperm diversification, such as pollination shifts (Kay & Sargent, 2009).
Different populations of plants occurring along elevational gradients can exhibit floral phenotypic variation mediated by pollinator selection (Galen, 1996;Maad et al., 2013;Nattero et al., 2011;Zhao & Wang, 2015). Studies of multiple bee-pollinated species found an increase in flower size with elevation (Galen, 1996;Maad et al., 2013;Malo & Baonza, 2002). Similarly, bee community composition can change with elevation, resulting in an increase in the mean body size in the community (Hoiss et al., 2012;Malo & Baonza, 2002). Research on the association between floral traits and bees' pollinator assemblages concluded that selection drives the increase of flower size along elevational gradients to improve the fit of bigger pollinators with the floral reproductive structures (Galen, 1996;Maad et al., 2013). Changes in pollinator assemblages can result in rapid floral trait divergence over a single generation (Galen, 1996) or even from one year to the next one (Schemske & Horvitz, 1989). To date, however, few studies have focused on floral traits and pollinator assemblages variation within species along elevational gradients and even fewer studies have been conducted in tropical areas (but see Dellinger et al., 2021 for a study focused on bee-to-hummingbird pollination shifts; and Klomberg et al., 2022 andNattero et al., 2011 for studies on geographic floral trait variation without pollination shifts). It is imperative to devote efforts to existing research programs on tropical species, such as species in
Bee diversity and abundance decrease with elevation in both temperate and tropical areas, (Arroyo et al., 1982;Hoiss et al., 2012).
In contrast, vertebrate pollinators, such as hummingbirds, can inhabit a wider elevational range than bees and perform well as pollinators. For instance, in rainy conditions at high elevations in Mexico, bird-pollinated plants were more effectively pollinated than closely related bee-pollinated plants (Cruden, 1972). The abundance and diversity of hummingbirds are higher in tropical versus temperate areas (Greenwalt, 1960). Thus, it is more likely that in the tropics, pollinator assemblages at high-elevation sites might include the presence of hummingbirds as pollinators, in addition to larger bee pollinators. There is however little empirical evidence regarding environmental factors, such as those associated with elevation (Thomson & Wilson, 2008), driving specific changes in pollinator assemblage from bees to hummingbirds (but see Dellinger et al., 2021), or from bees to a mixed pollinator assemblage.
Costus (Costaceae), a species-rich genus of herbaceous plants, has several Neotropical species pollinated either by orchid bees (Apidae: Euglossini) or hermit hummingbirds (Phaethornithinae, Kay & Schemske, 2003). Bee-and hummingbird-pollinated species have different suites of floral traits (Kay & Grossenbacher, 2022). Beepollinated species generally have long, wide, and pale flowers, and a large petaloid labellum colored with two thick yellow lines that presumably act as nectar guides. In contrast, hummingbird-pollinated species have short, narrow, and brightly colored flowers with a reduced labellum (Kay & Grossenbacher, 2022). Flowers of Costus guanaiensis var. tarmicus seem to differ in visual floral traits among sites distributed along an elevational gradient from the Amazonian basin to the Andes mountains in Peru (up to 2000 meters above sea level [m a.s.l.]). These observations and the broad elevational distribution of C. guanaiensis var. tarmicus offer a natural setup to study floral divergence across a tropical elevational gradient. We addressed the following questions: (1) Do floral traits vary among sites along an elevational gradient? (2) Do pollinator assemblages vary among sites along an elevational gradient? and (3) Is the pollinator assemblage variation associated with floral trait variation? To answer these questions, we conducted an analysis of C. guanaiensis var. tarmicus floral traits (size, color, and reward) and its pollinator assemblages at four sites along a 1000 m elevational gradient in Peru.
We had two hypotheses about how floral traits and pollinator assemblages of C. guanaiensis var. tarmicus might vary with elevation. First, if C. guanaiensis var. tarmicus is adapted to pollination by the local bee fauna (local bee adaptation hypothesis), we expected that the floral traits and pollinator assemblages would covary among sites. Given that there may be fewer bees in the highest site, and thus, a lower bee visitation rate compared to lower sites, attraction traits might be exaggerated there; for instance, the yellow lines (nectar guides) in the labellum would cover a wider area and the nectar sugar concentration would be greater at the high site compared to lower sites. Pollinator assemblages at all sites would comprise orchid bees, but the highest elevation site might receive visits from bigger bee species than those at lower elevation sites (Bishop & Armbruster, 1999;Galen, 1996;Maad et al., 2013). If so, we expect that the highest elevation site would have bigger flowers than lower elevation sites.
Alternatively, if in the high elevation site, C. guanaiensis var. tarmicus is not only adapted to the local bee fauna but also to hummingbird pollination (pollinator shift hypothesis), we expect that the floral traits would vary with elevation, as we explained in the previous hypothesis, but that the floral traits at the highest site would include traits that deter bees and attract or fit hummingbirds. Thus, we would not detect a flower and bee morphology correlation along the elevational gradient. For example, for a better hummingbird-flower fit, the flowers at the highest site would be the shortest to allow hummingbirds to reach the nectar and the narrowest to improve pollen transfer. To deter the bees, the flowers at the highest site would have thinner or absent yellow lines (nectar guides) on the labellum.
In addition, the nectar sugar concentration would be the lowest at the highest site, matching similar levels to known hummingbird flower sources (Baker, 1975;Bolten & Feinsinger, 1978). A dilute nectar would deter bees by reducing their sugar intake rate (Cnaani et al., 2006;Harder, 1986;Heinrich, 1975). Pollinator assemblages at the highest site would comprise orchid bees and hummingbirds, but with the lowest bee visitation rate and the highest hummingbird visitation rate compared to lower sites.

| Study system and sites
Costaceae is a family of monocots native to tropical climates of Central America, South America, Asia, and Africa (Maas, 1972). The studied species belongs to the Neotropical Costus clade with species distributed from Mexico to Brazil, encompassing the Andes mountainous regions, with species present in low-to mid-elevation (0-2000 m a.s.l.; Maas, 1972;Vargas et al., 2020). Costus typically produce a single nectar-rich flower per plant per day, and pollinators are thought to travel on traplines between widely spaced plants.
Bee pollination is ancestral in the Neotropical clade, but there have been at least 11 independent shifts to hummingbird pollination (Kay & Grossenbacher, 2022;Vargas et al., 2020), motivating the hypothesis that geographic divergence in floral traits might be associated with adaptation to hummingbirds.
Costus guanaiensis is a polyphyletic species with multiple named varieties (Maas, 1972;Vargas et al., 2020). Here we focus solely on the taxon C. guanaiensis var. tarmicus. Costus guanaiensis var. tarmicus ranges from 250 to 2000 m a.s.l. in the Peruvian Amazon and Andes (www.tropi cos.org). There are no previous studies on the pollination biology of this specific taxon, although there are records of orchid bees pollinating C. guanaiensis var. macrostobilus (Schemske, 1981;Sytsma & Pippen, 1985). The common name of these plants in Central and South America is "caña agria," due to their resemblance with sugar cane plants (Maas, 1972). (high elevation site, hereafter Yanachaga_H). During the 2019-2020 sampling season, we worked at sites Yanesha_L, Sanmatias_M, and Yanachaga_H but concluded that the logistics of working at Yanesha_L were too complicated to go back to the following sampling season.
Then, we decided to sample the Iscozacin_L site during 2020-2021 because we found a population of C. guanaiensis var. tarmicus there (not reported in the herbaria) more easily accessible than Yanesha_L. For the analysis, we decided to keep all four sites due to habitat conservation differences between the two low sites. All the sites are in the Department of Pasco on the eastern side of the Andes and western Amazonia in central Peru ( Figure 1).

| Floral phenotype analysis
We measured flower size, color corresponding to the nectar guides area, and nectar sugar concentration to quantify variation in floral traits among sites. We collected 9-11 flowers at each site to analyze flower size. In the morning, we searched for fresh flowers, removed them from the plants, and transported them in 50 mL thick plastic tubes to the research station/hotel, where we had a photography station set up. If any flower withered or damaged during the transportation, we did not use it for photos and waited to collect new flowers the next day. We took pictures of the flowers standing in molding clay to achieve a straight axis parallel to the camera. We took pictures of the frontal, lateral, and internal faces of each flower next to a ruler. From the pictures, we measured (in mm) 12 different floral traits using ImageJ (Rasband, 1997(Rasband, -2018. From the frontal view, we measured labellum length, labellum width, and corolla width. From the lateral view, we measured labellum length, corolla length, and corolla width. From the internal view, we measured petaloid stamen length (measured above the stigma), internal and horizontal labellum width, internal and diagonal labellum width, the distance between the anthers and the labellum, the distance between the anthers and the corolla base, and corolla tube length (see Appendix 1 for floral traits reference).
We performed a Principal Components Analysis (PCA) on floral traits measurements using the prcomp function in R (R Core Team, 2020). To assess if the floral traits differ among the sites, we applied a PERMANOVA followed by multiple pairwise comparison tests using the functions adonis (Thioulouse et al., 2018) and pairwise.adonis (Martinez Arbizu, 2020) from the package "vegan," respectively. Also, we assessed how the flowers vary specifically in length and width. To compare the floral length, we applied ANOVA tests to the frontal and lateral measurements of labellum length and the lateral measurement of corolla length. For the floral width, we applied ANOVA tests to the frontal and lateral measurements of the corolla width. Following ANOVAs, we analyzed pairwise comparisons using the post hoc test function TukeyHSD.
To quantify the colors corresponding to the nectar guides area, we took pictures of five flowers from Iscozacin_L, Sanmatias_M, and Yanachaga_H sites in 2021. We could not take pictures of Yanesha_L flowers due to the site's inaccessibility during our second sampling season. We took pictures of the top part of the floral labellum. We placed the flower next to a Spyderckr color card that served as the gray standard. We used a Sony Alpha R7 full spectrum converted camera, a UV transmission lens model Nikon EL-Nikkor 80mm f5.6, a visible-spectrum-pass filter (UV/IR Cut Hot Mirror 39 mm lens filter, Kolari Vision, US), and a UV-pass filter (UV Bandpass 39 mm lens filter, Kolari Vision, US). We took pictures under natural daylight between 1000 and 1200 h, avoiding moments when clouds were covering direct sunlight. We took several pictures with each filter varying the exposure levels to select the best shot for the analysis. We assessed the nectar guides area on the labellum using the MicaToolbox plugin (Troscianko & Stevens, 2015) in ImageJ.
First, we generated a multispectral image using one flower picture taken with the visible-spectrum-pass filter and one taken with the UV-pass filter. Then, we converted the multispectral image to a honeybee cone catch image (using the cone catch model available in the MicaToolbox software). Next, we acquired a presentation image for human vision from which we can distinguish colors seen by bees but not humans. These colors include the UV reflectance ones that we assumed corresponded to the flower's nectar guides.
In the image, they appeared in dark pink-purple colors. We converted this image to RGB colors and saved it as a jpeg image. Then, we converted the jpeg image to an 8-bit color type with 256 colors (maximum number available). By doing so, we aimed to better define the different color areas, especially those with dark pink-purple colors. Next, we converted the generated image to an 8-bit grayscale type to measure the dark pink-purple nectar guides area using the threshold settings. Then, we adjusted the threshold levels to identify the pixels corresponding to the nectar guides area. Finally, we used the measure command to get the nectar guides fraction area on the labellum (as the percentage area occupied on the labellum, see Appendix 2 for examples of the image processing). We compared the nectar guides fraction area among sites with a Kruskal-Wallis test, followed by a post hoc Dunn's test (due to the small sample size) for pairwise comparisons.
To quantify nectar sugar concentration, we collected nectar from 10 to 13 fresh flowers per site, between 0700 and 1030 h, using a capillary tube. The flowers we used were bagged the previous afternoon to avoid floral visitors taking the nectar before us. We placed a drop of the nectar in a refractometer (Eclipse refractometer, Bellingham + Stanley, UK) to record the sugar concentration in degrees Brix. We compared nectar sugar concentration among sites with an ANOVA test followed by Tukey's post hoc comparisons.

| Pollinator assemblage analysis
We recorded pollinators visiting the flowers of C. guanaiensis var.
tarmicus. We placed cameras with motion detection systems next to the inflorescence for 1-10 h per day. We used Canon Powershot To facilitate the identification of orchid bees in our videos, we caught orchid bees with entomological nets between 0900 and 1600 h for 3 days at each site. We always caught the orchid bees after we finished recording their pollination activity on the Costus flowers so that we did not interfere with their activity. We used pure chemical attractants eugenol, cineol, methyl salicylate, and benzyl acetate to attract them. We applied 0.10 mL of attractant to a cotton ball every half hour, and we used one cotton ball per attractant. We suspended the cotton balls at 1.50 m from the ground, 2 m apart among attractants. The captured individuals were euthanized and preserved with 70% ethanol for later identification in the Entomology Laboratory of the Universidad Nacional San Antonio Abad del Cusco, Peru, using a stereoscope (NOVEL NSZ-608T). For the determination of the specimens, we used the terminology and taxonomic keys proposed by Dressler (1984), Bonilla-Gomez and Nates-Parra (1992), and Michener (2000). We tried to determine the orchid bee in our videos to species level when possible. Moreover, we categorized bees into three groups based on their size, small (Euglossa spp), medium (Eulaema cf mocsaryi, Eulaema cf polychroma, Eulaema sp.), and big (Eulaema cf bombiformis, Eufriesea cf ornata), to estimate visitation rate by bee size. We decided to exclude the species Aglae caerulea from the big-size bee group because it belongs to a different genus which is known to have ecological differences with species of the genus Eulaema (Roubik & Hanson, 2004) and we only recorded three visits of A. caerulea at Yanachaga_H site.
We also constructed a standardized variable called pollinator morphotype abundance to compare pollinator assemblages among sites.
This was calculated by dividing the visits per flower of each pollinator by the total number of pollinators visits at each site.
For the statistical analysis, we first compared the visitation rate of each pollinator group (bees separately from hummingbirds) among sites with a non-parametric Kruskal-Wallis test, and a pairwise Wilcoxon post hoc test. Second, we compared the visitation rate of different orchid bee sizes among sites with a Kruskal-Wallis test and a pairwise Wilcoxon post hoc test. The visitation rate of big-sized orchid bees was mostly zero in mid and low sites, then, for the Wilcoxon test, we indicated the alternative hypothesis of "less" for the group that had all or almost all zero values. Third, we compared the bee visitation rate by size within each site by applying the Kruskal-Wallis and the post hoc pairwise Wilcoxon tests. Fourth, we used the pollinator morphotype abundance to perform a Non-metric Multidimensional Scaling (NMDS) using the metaMDS function in the package "vegan" (Oksanen et al., 2020). We applied a PERMANOVA and a multiple pairwise comparison test, to assess if pollinator assemblages differ among sites. Also, we ran a Simper analysis using the function simper from the package "vegan" to identify the pollinator morphotypes from each site that contributed to the dissimilarity with the other sites. Finally, from the bee sampling, we estimated orchid bee relative abundance in the different sites to evaluate if the orchid bee abundance changes along the elevational gradient. We applied a Kruskal-Wallis test and the post hoc Dunn's test (due to small and unequal sample size among the sites).
In addition, we carried out a survey of hummingbirds in the sites Iscozacin_L, Sanmatias_M, and Yanachaga_H during the 2020-2021 season. We conducted four walk transects of 1 km in each site and recorded the hummingbird species by observations (using binoculars) or songs. We looked up the bill length of the recorded hummingbird species in the literature (Hilty, 2002;Meyer et al., 1970;Schulenberg et al., 2010) to assess if they could reach the nectar of the flowers of C. guanaiensis var. tarmicus in the different sites. We considered hummingbird species that conducted a legitime visit to the flowers as having long enough bill length to reach the nectar of the flowers. We aimed to identify if there were other species with similar bill length in the sampled sites.

| Floral phenotypes and pollinators size association
We determined whether floral size traits correlated with orchid bee body traits to evaluate if the different floral phenotypes were associated with the pollinator sizes. We chose the bee species with the highest visitation rate within each site to take the bee body measurements, assuming that the species most likely effecting the highest selection pressure on the floral traits would be the most frequent pollinator (Stebbins, 1970). We selected the species Euglossa imperialis in Iscozacin_L, Euglossa intersecta in Yanesha_L, Eulaema mocsaryi in Sanmatias_M, and Eulaema bombiformis in Yanachaga_H. We measured the bee thorax width and height, and the bee tongue (proboscis) length from 5 to 10 individuals of each orchid bee species.
Using the mean of each variable, we applied Pearson's correlation between the following bee and flower traits, bee thorax width and the frontal corolla width, bee thorax height and the lateral corolla width, bee tongue length and corolla tube length (or the functional floral length that interacts with the bee's tongue), and bee tongue length and the distance between the anthers and the corolla base.
We also conducted a bootstrap resampling from our empirical bee body measurements with replacement to pair each individual plant measurement with an individual bee measurement per site, thus, incorporating variation in each trait. We repeated the resampling for TA B L E 1 Pollinator counts and visitation rate per hour to inflorescences of four populations of Costus guanaiensis var. tarmicus along an elevational gradient.

| Floral phenotypes vary among sites
We found that floral traits varied among sites. Flower size and nectar sugar concentration differed among two or more sites. The We found variation in nectar guides fraction area among sites (Kruskal-Wallis, χ 2 = 11.02, df = 2, p = .004). Dunn's test revealed that Yanachaga_H was significantly different than Sanmatias_M (p = .002) and Iscozacin_L (p = .05). But the latter two did not differ between them. Thus, the nectar guides fraction area was the smallest in Yanachaga_H (x̄ = 0.37%, SD ± 0.83, Figure 3a). Similarly, we found variation in nectar sugar concentration among sites (ANOVA, F = 6.47, p < .001). The Tukey's post hoc test showed that Yanachaga_H was significantly different from the other three sites (p < .05). Thus, Yanachaga_H had the lowest nectar sugar concentration (x̄ = 38.9°B, SD ± 1.4), whereas the other three sites had similar levels of nectar sugar concentration (Figure 3b).

| Pollinator assemblages vary among sites
We video-recorded a total of 1192 h of pollinator activity on the flowers of C. guanaiensis var. tarmicus. The flowers were visited by bees and hummingbirds ( However, the hummingbird visitation rate is only approximately 7% of the bee visitation rate in Yanachaga_H (Table 1).
We compared the visitation rates of small, medium, and big-sized orchid bees among sites (Figure 4a). We found a significant differ-

| Floral phenotype and pollinator size association
We tested if floral traits and bee traits covaried using Pearson's correlation test. We found a non-significant correlation between our site-averaged variables, as expected with four sites. The correlation between the bee body size traits and the floral width traits showed positive correlation coefficients higher than .75 ( Figure 5a,b), whereas the correlation between the tongue length and the corolla tube length, and the tongue length and the distance between the anthers and the corolla base showed a negative correlation coefficient lower than .4 (Figure 5c). After bootstrap resampling and repeated correlation, we found a positive significant correlation between bee thorax width and the frontal corolla width (1000 out of 1000 repetitions with a p < .05, .34 < coef < .5), and between the bee thorax height and the lateral corolla width (1000 out of 1000 repetitions with a p < .05, .42 < coef < .57). We did not find a significant correlation between the bee tongue length and the corolla tube length (61 out of 1000 repetitions with a p < .05), nor between the bee tongue length and the distance between the anthers and the corolla base (0 repetitions out of 1000 with a p < .05).

| DISCUSS ION
Our results showed clear variation in floral traits and pollinator assemblages in C. guanaiensis var. tarmicus across a steep elevational gradient from the Amazon to the foothills of the Andes in Peru.
We also found an association between floral traits and bee body traits contributing to the mechanical flower-pollinator fit. The substantial divergence in both floral traits and pollinator assemblages suggests the establishment of the necessary preconditions for pollinator-driven divergence. Our findings highlight the TA B L E 2 p-Values from the pairwise Wilcoxon tests comparing orchid bee visitation rate by size within each sampling site in our elevational gradient.

F I G U R E 5
The association between the bee body traits of the species with the highest visitation rate per site and the flower size traits. (a) A positive correlation between the bee thorax width and the frontal corolla width. (b) A positive correlation between the bee thorax height and the lateral corolla width. (c) No correlation between the bee tongue length and the corolla tube length, and between the bee tongue length and the distance between the anthers and the corolla base.
importance of assessing floral traits and pollinator assemblages' geographic variation in the Neotropics by allowing us to learn from a new plant study system and novel composition of pollinator assemblages that are not found in temperate areas (orchid bees and hermit hummingbirds). Below we examine our results in light of two hypotheses for how variation is structured: that it could reflect local adaptation to different assemblages of the same pollinator functional group (bees) or that it could reflect the early stages of an evolutionary shift to a new pollinator functional group (hummingbirds).

| Is there evidence for the local bee adaptation hypothesis?
The local bee adaptation hypothesis predicts that the pollinator as-  Galen, 1996; and Campanula rotundifolia in the Norwegian mountains by Maad et al., 2013). Previous studies found an association between the flower length and the pollinator mouth part length (Maad et al., 2013;Nagano et al., 2014), which we did not find in our study. Perhaps the flower length is not a trait upon which orchid bees' tongue length exerts a selective pressure, given that the bees crawl inside the flowers of Costus to feed on the nectar. Studies about the evolution of nectar flowers for orchid bees suggested that long flowers for orchid bees might have evolved via competition among sympatric species for attracting traplining pollinators that would include them in their feeding routes (Borrell, 2005;Garrison & Gass, 1999;Rathcke, 1992), and not via directional selection exerted by specialized pollinators (Darwin, 1862), such as those with extreme trait adaptations. In addition, the local bee adaptation hypothesis stated that if there is a low bee visitation rate in the highest site, we should see an exaggeration of bee attraction traits there. Our results did not support this part of the hypothesis; we neither found a lower bee visitation rate nor exaggerated bee attraction traits in the highest site. This finding contrasts with the study of Dellinger et al. (2021), who found a decrease in bee visitation rate with elevation in Merianieae species, but their high elevation sites were located beyond 1800 m a.s.l.
Altogether, we found evidence that supports the adaptation of C. guanaiensis var. tarmicus flowers to the local bee fauna along the studied elevational gradient.

| Is there evidence for the hummingbird pollination shift hypothesis?
The pollinator shift hypothesis predicts that the pollinator assem-  (Biesmeijer et al., 2006;Maglianesi et al., 2015;Ornelas et al., 2007;Smith et al., 1995). Regarding hummingbird fitting floral traits, the short corolla in the flowers of the highest site should make those flowers accessible to the hummingbirds. In addition, the level of sugar concentration in the nectar of the highest site (x̄ = 38.9° Brix, SD ± 1.4) might work as a hummingbird attracting floral trait. This is because the sugar level in the highest site is higher than the sugar level found in many hummingbird-visited flowers (Baker, 1975;Bolten & Feinsinger, 1978;Chalcoff et al., 2006;McDade & Weeks, 2004;Rodríguez-Flores & Stiles, 2005), including hummingbird-pollinated Costus species (Rodríguez-Flores & Stiles, 2005;Sytsma & Pippen, 1985). We initially thought that the level of sugar concentration in the nectar found in the highest site might be related to bee deterrence since it is lower than the one found in the lower sites. However, the sugar level found in the nectar of the highest site is within the range of sugar level of most flowers foraged by euglossine bees (30%-40%, Roubik et al., 1995). Regarding the availability of other hummingbird floral resources, unfortunately, we did not account for other hummingbird floral resources in our studied sites to analyze if there are lesser available resources in the highest site compared to lower sites. Maglianesi et al. (2015) showed that if there is a decrease in hummingbird floral resources, hummingbirds can show a less specialized diet feeding on a larger number of plant species. Moreover, Stiles (1978) reported hermit hummingbirds feeding on a bee-pollinated Costus species (C. malortieanus) in Costa Rica at the end of the rainy season (November-December) which matches the low flowering point of hummingbird foodplants during the year.
Another finding supporting the pollination shift hypothesis is the reduced nectar guides fraction area in the highest site compared to lower sites. Schemske and Bradshaw (1999) found that the removal of nectar guides on the Mimulus flowers lowered the bee visitation rate. In our study, we did not find a decrease in bee visitation rate, perhaps the nectar guides area must be completely absent to influence bee visitation. Overall, we found some evidence supporting the pollinator shift hypothesis. Although the association between the corolla size and the bee thorax size of the most frequent pollinator could contradict this hypothesis, the two scenarios might not be mutually exclusive.

| Floral traits and pollinator assemblages differed between the lowest sites
Flower length differed between the two lowest sites, Iscozacin_L has shorter flowers than Yanesha_L. Moreover, we also found differences in pollinator assemblages and pollinator visitation rates between them. These findings suggest that other factors beyond the elevation can promote variation in floral traits and pollinator assemblages.
Furthermore, the Yanesha_L bee's visitation rate was the highest among all the sites. This result might be driven by differences in the level of habitat degradation among sites. Yanesha_L is in the middle of protected areas, whereas Iscozacin_L is outside of a protected area surrounding a small town. The other sites, Sanmatias_M and Yanachaga_H, are on the border of protected areas with roads crossing these sites (Figure 1). Previous studies showed that orchid bee abundance and richness are higher in well-preserved areas compared to disturbed areas (Storck-Tonon & Peres, 2017). Our findings contribute evidence to the importance of habitat conservation for orchid bees and their role as pollinators, considering that, in average, the bee visitation rate inside a protected area was 7 times higher compared to unprotected areas. This striking difference may also affect the reproductive success of our native studied plant species C. guanaiensis var. tarmicus.

| Correlation or causation? Caveats about the association of floral traits and pollinator assemblages
Some caveats about correlation versus causation for the association of floral traits and pollinator assemblages merit discussion. First, we do not know the mechanism behind the correlation we found between the corolla size and the bee thorax size. The flower size may adapt to the available pollinators, or the flower size may vary due to pollination-unrelated factors and the pollinators preferentially feed on the right size of flowers for them (Nagano et al., 2014). If the latter mechanism was true, we might see

| CON CLUS IONS
We found substantial divergence in floral traits of a Neotropical Costus species and its pollinator assemblage across a steep elevational gradient spanning the Amazonian lowlands to the eastern foothills of the Central Andes, establishing the necessary preconditions for pollinator-driven divergence. Taking together the results from floral traits and pollinator assemblages variation suggest that the populations of C. guanaiensis var. tarmicus are adapted to the local bee fauna along the studied elevational gradient, but we cannot rule out the possibility of the beginning of a bee-to-hummingbird pollination shift in the highest studied site.

ACK N OWLED G M ENTS
We thank our funding sources, including the Faculty for the Future

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data analyzed in the manuscript will be made accessible through the public Dryad repository. https://datad ryad.org/stash/ share/ _Y74fA aqB4m NN4rW ugtTt osIUC Y4aaB Q w50A R1tf6Os.